DE between covid and control by cell type

Started: 2020-09-25
Last edited: 2020-10-01 13:31:50

library(tidyverse)

# single cell
library(Seurat)

# plotting
library(patchwork)
library(ggthemes)
library(ggrepel)
library(ggforce) # for sina plots
library(cowplot)
library(ggtext)
library(ungeviz) # for geom_hpline

# go enrichment
library(clusterProfiler)
library(org.Hs.eg.db)
library(GO.db)

1 Background

Our goal here is to perform differential expression analysis between covid and control samples separately for each annotated cell type in the single cell RNA-seq data.

This analysis is the same as 05_de-by-celltype_logfc0.25.Rmd, except for the log fold change cutoff for DE analysis.

2 Data

We will use the Seurat object in which we have saved the cell annotations.

seur <- readRDS("../results/02_annotation/seurat-object_annotated.rds")

3 Differentialy expressed genes

We will use the “merged” annotations as our celltypes for which to perform DE analysis. For identifying differentially expressed genes, we will use Seurat::FindMarkers(). For this we first have to create a new metadata column with celltype_covid and celltype_control values, and then we run DE between these two idents. This function by dafault uses the Wilcoxon Rank Sum test, which we won’t change. By default this function outputs all genes for which abs(logFC) is at least 0.25, without filtering by p value. We will not use these defaults. Instead:

1. Fold change cutoff: \(|log(Fold Change)| > 0.4\)
   
2. \(p_{adj} < 0.05\)

Seurat::FindMarkers() uses natural log, not log2, and uses bonferroni adjustment to correct for multiple testing (see https://github.com/satijalab/seurat/issues/741). \(|log(FoldChange)| > 0.4\) corresponds to exp(0.4) = 1.49 or at least 50% change in expression.

The output of this for each cell type is stored as .csv in results/04_de-genes-by-celltype/logfc_0.40/de/files.

# set up
seur$celltype_covid <- paste(seur$annotation_merged, seur$covid, sep = "_")
Idents(seur) <- seur$celltype_covid

# de
de.genes <- list()

for (i in unique(seur$annotation_merged)) {
  de.genes[[i]] <-  FindMarkers(seur, 
                                ident.1 = paste0(i, "_covid"), 
                                ident.2 = paste0(i, "_cntrl"),
                                assay = "SCT", verbose = FALSE) %>% 
    dplyr::filter(p_val_adj < 0.05 & abs(avg_logFC) > 0.4) %>% 
    dplyr::rename(pct.covid = pct.1,
                  pct.cntrl = pct.2)
}

# write DE genes to files
for(i in unique(seur$annotation_merged)) {
  write.csv(
    de.genes[[i]],
    file = paste0("../results/04_de-genes-by-celltype/logfc_0.40/de/files/de-genes_", i, ".csv" ),
    row.names = TRUE)
}

3.1 Number of DE genes by celltype

de.number <- lapply(de.genes, nrow) %>% as.data.frame() %>% t() %>% as.data.frame()
colnames(de.number) <- "number"
de.number$celltype <- rownames(de.number)

de.number <- de.number[order(de.number$number, decreasing = TRUE), ]
de.number$celltype <- factor(de.number$celltype, levels = de.number$celltype)

p <- ggplot(data = de.number, aes(x = celltype, y = number, label = number)) +
  geom_bar(stat = "identity", fill = "grey50", width = 0.8) +
  geom_text(size = 1.8, vjust = -0.25) +
  ggtitle("Number of DE genes") +
  coord_cartesian(clip = "off") +
  theme_minimal() +
  theme(plot.title = element_text(size = 8),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(size = 0.2),
        legend.key.size = unit(0.75, "lines"),
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6),
        axis.text.y = element_text(size  = 6),
        axis.title = element_text(size = 7))

ggsave(p, filename = "../results/04_de-genes-by-celltype/logfc_0.40/de/plots/number-of-de-genes.pdf", width = 3.5, height = 3, units = "in")

print(p)

3.2 Mitochondrial UMI

For some cell types a lot of mitochondrial genes are DE between covid and control. It’s possible that those genes are really differentially expressed, but it’s more likely that there is a difference between covid and control samples in terms of quality of cell-prep. We can look at percent mitochondrial UMI in all cell types by covid status to see if this is true.

ggplot(data = seur@meta.data, aes(x = covid, y = percent.mt)) +
  geom_sina(aes(color = covid), size = 1, alpha = 1) +
  scale_color_tableau() +
  facet_grid(.~ annotation_merged) +
  labs(caption = "% mitochondrial reads in all cell types") +
  theme_minimal() +
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        panel.grid.major = element_line(size = 0.2),
        panel.grid.minor = element_blank(),
        strip.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5)) 

For most cell types there appears to be a difference in the distribution of percent mitochondrial reads between covid and control samples. Covid samples generally have lower % mitochondrial reads. Next, we can plot the absolute value of the difference in percent mitochondrial reads between covid and control against the number of DE genes that are mitochondrial. If the mitochondrial genes in DE analysis are driven by differences in sample quality, we should expect a positive relationship.

mt <- data.frame(
  celltype = sort(names(de.genes)),
  mt.diff = NA,
  mt.de = NA
)

for(i in mt$celltype) {
  mt$mt.diff[mt$celltype == i] <- abs(
    median(seur$percent.mt[seur$annotation_merged == i & seur$covid == "covid"]) -
      median(seur$percent.mt[seur$annotation_merged == i & seur$covid == "cntrl"])
  )
  
  mt$mt.de[mt$celltype == i] <- grep("MT-", rownames(de.genes[[i]])) %>% length()
}

ggplot(data = mt, aes(x = mt.diff, y = mt.de, label = celltype)) +
  geom_point(shape = 21, color = "black", fill = "grey40", size = 3) +
  coord_cartesian(xlim = c(0, NA), ylim = c(0, NA)) +
  geom_text_repel(size = 3, color = "grey40") +
  labs(x = "\nabs(median(% mt in covid) - median(% mt in control))",
       y = "number of DE mitochondrial genes") +
  theme_minimal() +
  theme(panel.grid.major = element_line(size = 0.25),
        panel.grid.minor = element_blank())

There is a very mild positive relationship, more so for villous cell types. Three decidual cell types—dec.Mono_1, dec.Mono_2, and dec.Tcell_2—have high number of MT DE genes despite low difference in percent MT reads.

We can ignore the mitochondrial genes since the number of genes is quite small compared to the total number of DE genes.

4 Interferome

Interferome (Rusinova et al. 2012) is a manually curated database of interferon-responsive genes, available at http://www.interferome.org/interferome/home.jspx. The database essentially contains a list of all genes whose expression is affected by interferons, as inferred from various datasets. These datasets include in vitro and in vivo studies from human and mouse. The online interface includes a search function which allows users to input a list of genes, which are searched against the database. The output includes various things including: (1) list of genes from the input that are found in the inferferome database, i.e. there is evidence for these genes to be modulated by interferon, (2) a detailed table with links pointing to the dataset from which the evidence was gathered for each gene, (3) gene ontology enrichment, (4) TF binding site analysis, etc.

I ran the list of DE genes for each of the 21 celltypes through Interferome (version 2.01, 2020-09-20) and saved the text output for (1) from above. I used the following search parameters (all defaults except for narrowing species to Homo sapiens)

Interferome Type:   Any
Interferome SubType:    Any
Treatment Concentration:    Any
Treatment Time: Any
Vivo/Vitro: Any
Species:    Homo sapiens
System: Any
Organ:  Any
Cell:   Any
Cell Line:  Any
Normal/Abnormal:    Any
Fold Change Up: 2.0
Fold Change Down:   2.0

4.1 Fraction of DE in interferome database

We want to calculate for each celltype, the fraction of DE genes (between covid and control) which are known to be interferon-responsive. For this, below we will read the Interferome output text files and calculate that fraction.

# read interferome output
ifome <- list()

for(i in names(de.genes)) {
  ifome[[i]] <- read.table(
    paste0(
      "../results/04_de-genes-by-celltype/logfc_0.40/interferome/files/interferome_", 
      i, ".txt"), 
    header = TRUE, sep = "\t", quote = "", skip = 19, 
    strip.white = TRUE, stringsAsFactors = FALSE)
}

# create a dataframe of number of genes in interferome for each cell type
ifome.n <- data.frame(
  celltype = names(de.genes),
  n_de = NA,
  n_interferome.yes = NA
)

# calculate number of genes
for (i in names(de.genes)) {
  
  ifome.n$n_de[ifome.n$celltype == i] <- nrow(de.genes[[i]])
  
  ifome.n$n_interferome.yes[ifome.n$celltype == i] <- 
    (rownames(de.genes[[i]]) %in% ifome[[i]]$Gene.Name) %>% 
    sum()

}

ifome.n$n_interferome.no <- ifome.n$n_de - ifome.n$n_interferome.yes

# calculate percentages
ifome.n$pct_interferome.yes <- (ifome.n$n_interferome.yes/ifome.n$n_de) * 100
ifome.n$pct_interferome.no <- 100 - ifome.n$pct_interferome.yes

# prepare for plotting
# melt
ifome.n.long <- pivot_longer(data = ifome.n,
                             cols = 3:6,
                             names_to = c(".value", "interferome"), 
                             names_pattern = "(.+)_(.+)")
ifome.n.long$interferome <- gsub("interferome.", "", ifome.n.long$interferome)

ifome.n.long$celltype <- factor(ifome.n.long$celltype, 
                                levels = de.number$celltype)

# plot iterferome output

theme2 <-   theme_classic() +
  theme(
    legend.key.size = unit(0.5, "lines"),
    legend.title = element_text(size = 6),
    legend.text = element_text(size = 5.5),
    axis.text.y = element_text(angle = 0, hjust = 1, vjust = 0.5, size = 6, 
                               color = "black"),
    axis.text.x = element_text(size = 6, color = "black"),
    axis.title = element_text(size = 7),
    axis.title.y = element_blank(),
    axis.line = element_line(size = 0.25),
    axis.ticks.y = element_line(size = 0.25),
    axis.ticks.x = element_line(size = 0.25))

p.n <- ggplot(data = ifome.n.long, 
              aes(y = celltype, x = n)) +
  geom_bar(position = "stack", stat = "identity", 
           aes(fill = interferome), width = 0.8) +
  scale_fill_manual(values = c("grey90", "#4e79a7"), 
                    breaks = c("no", "yes")) +
  scale_y_discrete(expand = c(0.05, 0)) +
  scale_x_continuous(expand = c(0.025, 0)) +
  geom_text(data = ifome.n.long[ifome.n.long$interferome == "yes", ], 
            aes(y = celltype, x = n, label = n),
            size = 1.8, vjust = 0.5, hjust = 0.0, nudge_x = 1) +
  labs(x = "number") +
  theme2 +
  theme(legend.position = c(0.90, 0.90))

p.pct <- ggplot(data = ifome.n.long, 
                aes(y = celltype, x = pct)) +
  geom_bar(position = "stack", stat = "identity", 
           aes(fill = interferome), width = 0.8) +
  scale_fill_manual(values = c("grey90", "#4e79a7"), 
                    breaks = c("no", "yes")) +
  scale_y_discrete(expand = c(0.05, 0)) +
  scale_x_continuous(expand = c(0.025, 0)) +
  geom_text(data = ifome.n.long[ifome.n.long$interferome == "yes", ], 
            aes(y = celltype, x = pct, label = round(pct, 0)),
            size = 1.8, vjust = 0.5, hjust = 0.0, nudge_x = 0.5) +
  labs(x = "percent") +
  theme2 +
  theme(legend.position = "none")


plot_row <- plot_grid(p.n, p.pct)

# now add the title
title <- ggdraw() +
  draw_label("DE genes in interferome database", x = 0, hjust = 0, size = 8) +
  theme(
    # add margin on the left of the drawing canvas,
    # so title is aligned with left edge of first plot
    plot.margin = margin(0, 0, 0, 7)
  )

p <- plot_grid(title, plot_row,
               ncol = 1,
               # rel_heights values control vertical title margins
               rel_heights = c(0.05, 1)
)

ggsave(p,
       filename = "../results/04_de-genes-by-celltype/logfc_0.40/interferome/plots/de-genes-interferome.pdf",
       width = 6, height = 3, units = "in")

print(p)

4.2 Test for enrichment

For most cell types, a large chunk of genes differentially expressed between covid and control are in fact known to be responsive to interferons according to Interferome. To find out whether this fraction is larger than what one would expect by chance, we can do a hypergeometric test for enrichment, similar to the one that we do for GO enrichment.

If \(N\) is the total number of genes used in DE analysis, \(B\) is the number of genes out of \(N\) that are in the interferome database, \(n\) is the number of DE genes for a cell type, and \(b\) is the number of genes genes out of \(n\) that are in the interferome database, then: \(enrichment = \frac{b}{n}/\frac{B}{N}\). And the hypergeometric distribution, we can calculate the p value (the cumulative probability of finding \(b\) or more genes in a sample of \(n\) DE genes, given \(B\) and \(N\)).

# total number of genes in the dataset
ifome.n$N <- nrow(seur@assays$SCT@data)

# total number of genes in the interferome dataset
ifome.n$B <- 12614

# enrichment
ifome.n$enrichment <- (ifome.n$n_interferome.yes/ ifome.n$n_de)/
  (ifome.n$B/ifome.n$N)

# p value
#phyper(q = 244-1, m = 12614, n = 23502-12614, k = 355, lower.tail = FALSE)

ifome.n$pval <- phyper(
  q = ifome.n$n_interferome.yes - 1,
  m = ifome.n$B,
  n = ifome.n$N - ifome.n$B,
  k = ifome.n$n_de,
  lower.tail = FALSE
)

# colors for plot
ifome.n$color <- ifelse(ifome.n$pval < 0.05, yes = "Black", no = "Grey")

# write to file
ifome.n %>% 
  dplyr::rename("n_total.genes" = "N", "n_interferome.genes" = "B") %>% 
  dplyr::select(-"color") %>% 
  write.csv(., "../results/04_de-genes-by-celltype/logfc_0.40/interferome/files/number-of-DE-genes-in-interferome.csv", 
            row.names = FALSE)

ifome.n[, c(1, 2, 3, 7, 8, 9, 10)] %>% knitr::kable()

celltype	n_de	n_interferome.yes	N	B	enrichment	pval
dec.DSC	355	244	23502	12614	1.280600	0.0000000
dec.SMC	276	204	23502	12614	1.377124	0.0000000
dec.APC	267	229	23502	12614	1.597998	0.0000000
vil.FB	161	133	23502	12614	1.539139	0.0000000
dec.Endo	263	191	23502	12614	1.353099	0.0000000
vil.Hofb	131	107	23502	12614	1.521824	0.0000000
dec.NK_1	348	261	23502	12614	1.397376	0.0000000
vil.EVT	61	40	23502	12614	1.221749	0.0401057
dec.FB	193	155	23502	12614	1.496327	0.0000000
vil.VCT	118	100	23502	12614	1.578956	0.0000000
dec.Tcell_1	349	253	23502	12614	1.350663	0.0000000
dec.Mono_1	176	125	23502	12614	1.323273	0.0000017
dec.Tcell_2	83	61	23502	12614	1.369316	0.0001628
dec.NK_2	276	206	23502	12614	1.390625	0.0000000
dec.Bcells	23	15	23502	12614	1.215110	0.1840788
vil.Ery	232	161	23502	12614	1.292974	0.0000007
dec.NK_3	282	201	23502	12614	1.328003	0.0000000
dec.Tcell_3	116	74	23502	12614	1.188573	0.0173755
dec.Mono_2	112	91	23502	12614	1.513824	0.0000000
dec.Gran	5	0	23502	12614	0.000000	1.0000000
vil.SCT	23	13	23502	12614	1.053095	0.4762411

# Plot fraction of DE genes in interferome, along with enrichment p values.

plot.dat <- ifome.n.long
plot.dat$celltype <- factor(
  plot.dat$celltype, 
  levels = ifome.n$celltype[order(ifome.n$pct_interferome.yes)]
  )
p <- ggplot(data = plot.dat, 
            aes(y = celltype, x = pct)) +
  geom_bar(position = "stack", 
           stat = "identity", 
           aes(fill = interferome), 
           width = 0.8) +
  scale_fill_manual(values = c("grey90", "#4e79a7"), 
                    breaks = c("no", "yes"),
                    name = "Interferome", 
                    labels = c("No", "Yes")) +
  scale_y_discrete(expand = c(0.05, 0)) +
  scale_x_continuous(expand = c(0.025, 0)) +
  coord_cartesian(clip = "off") +
  geom_text(data = plot.dat[plot.dat$interferome == "yes", ], 
            aes(y = celltype, x = pct, 
                label = paste0(n, "/", n_de)),
            size = 1.8, vjust = 0.5, hjust = 0.0, nudge_x = 0.5) +
  labs(x = "DE genes present in Interferome (%)") +
  geom_text(data = ifome.n[order(levels(plot.dat$celltype)), ], 
            aes(y = celltype, x = 100, 
                label = formatC(pval, format = "e", digits = 0)),
            hjust = 0, size = 1.75, nudge_x = 1,
            color = ifome.n$color[order(levels(plot.dat$celltype))]) +
  annotate(geom = "richtext", x  = 101, y = 22, size = 2,
           label = "enrichment *p* value", hjust = 0,
           fill = NA, label.color = NA, # remove background and outline
           label.padding = grid::unit(0, "pt")) + # remove padding
  theme2 +
  theme(legend.position = "right",
        legend.margin = margin(0, 0, 0, 3))

ggsave(p, filename = "../results/04_de-genes-by-celltype/logfc_0.40/interferome/plots/interferome_percent_enrichment-pval.pdf", width = 3.5, height = 3, units = "in")

print(p)

Even though the enrichment is modest, about 1.5 fold, for most cell types, it is highly unlikely by chance. The p values for most cell types are almost zero. Only three cell types have p values higher than 0.05.

5 GO enrichment

I also ran GO enrichment (only for the Biological Process category) for DE genes (only those that are up-regulated in covid) in all celltypes. The table outputs for all celltypes are saved as .csv in results/04_de-genes-by-celltype/logfc_0.40/go/files. For most cell types, the enriched GO categproes include those related to defense to viral infection, interferon response genes, misfolded protein response, etc. See below for the top 10 enriched GO categories for DE genes in each celltype.

# function for GO enrichment on upregulated genes for each cell type
enrichGO2 <- function(celltype, ...) {
  
  up <- de.genes[[celltype]] %>% 
    dplyr::filter(avg_logFC > 0)
  
  ids.bitr <- bitr(rownames(up), 
                   fromType = "SYMBOL", toType = c("ENTREZID", "SYMBOL"), 
                   OrgDb = "org.Hs.eg.db")
    
  features <- ids.bitr$ENTREZID

  ego <- enrichGO(
    gene          = features,    
    OrgDb         = org.Hs.eg.db,
    ont           = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff  = 0.01, # default 0.01
    qvalueCutoff  = 0.05, # default 0.05
    readable = TRUE,
    ...)
  
  dftoprint <- clusterProfiler::simplify(ego)@result
  rownames(dftoprint) <- NULL
  
  return(dftoprint)
}

# perform GO enrichment
de.go <- list()

for(i in names(de.genes)) {
  de.go[[i]] <- enrichGO2(celltype = i)
}

'select()' returned 1:1 mapping between keys and columns
2.09% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
1.93% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
0.81% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
2.46% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
1.27% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns
2.91% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns
2% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
1.89% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
4.17% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
2.38% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
3.57% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
3.38% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns
4.76% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
3.03% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
4.31% of input gene IDs are fail to map...'select()' returned 1:1 mapping between keys and columns
2.5% of input gene IDs are fail to map...--> No gene can be mapped....
--> Expected input gene ID: 51668,7455,1017,23493,7486,5266
--> return NULL...
Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘simplify’ for signature ‘"NULL"’

# write GO enrichments to files
for(i in unique(names(de.go))) {
  write.csv(
    de.go[[i]],
    file = paste0("../results/04_de-genes-by-celltype/logfc_0.40/go/files/enrichedGO_", i, ".csv" ),
    row.names = FALSE)
}

for(i in names(de.go)) {
  print(i)
  print(as.data.frame(de.go[[i]][1:10, 2]))
}

[1] "dec.DSC"

[1] "dec.SMC"

[1] "dec.APC"

[1] "vil.FB"

[1] "dec.Endo"

[1] "vil.Hofb"

[1] "dec.NK_1"

[1] "vil.EVT"

[1] "dec.FB"

[1] "vil.VCT"

[1] "dec.Tcell_1"

[1] "dec.Mono_1"

[1] "dec.Tcell_2"

[1] "dec.NK_2"

[1] "dec.Bcells"

[1] "vil.Ery"

[1] "dec.NK_3"

[1] "dec.Tcell_3"

[1] "dec.Mono_2"

6 Volcano plots

Below are volcano plots for all cell types: -log10(p value) vs logFC; top 30 genes with highest abs(logFC) are labeled.

# function for volcano plot
de_volcano <- function(celltype) {
  
  plot.dat <- de.genes[[celltype]]
  plot.dat$gene <- rownames(plot.dat)
  
  # top DE genes
  genes.to.label = de.genes[[celltype]] %>% 
    .[order(abs(de.genes[[celltype]]$avg_logFC), decreasing = TRUE), ] %>% 
    rownames() %>% .[1:30]
  
  plot.dat$label <- ifelse(test = plot.dat$gene %in% genes.to.label,
                           yes = plot.dat$gene,
                           no = NA)
  
  p <- ggplot(plot.dat, 
              aes(avg_logFC, -log10(p_val_adj), label = label)) + 
    geom_point(alpha = 0.4, size = 2) + 
    ggtitle(celltype) +
    geom_label_repel(size = 1.9, force = 1, color = "grey20") +
    theme_bw() +
    theme(panel.grid.major = element_line(size = 0.25),
          panel.grid.minor = element_blank(),
          axis.title = element_text(size = 7),
          axis.text = element_text(size = 7),
          plot.title = element_text(size = 8)) 
  
  return(p)
}

for(i in names(de.genes)) {
  p <- de_volcano(celltype = i)
  
  ggsave(
    p,
    filename = paste0(
      "../results/04_de-genes-by-celltype/logfc_0.40/de/plots/volcanoplot_", 
      i, ".pdf"),
    width = 5, height = 4, units = "in")
  
  print(p)
}

7 Labor signature

The control samples in our data come from C-sections. Which makes this a confounding factor in our analyses. In order to address this, we can compare the list of DE genes between covid and control samples with the list of DE genes between term in labor (TIL) and term not in labor (TNL). This can tell us if (and how much) part of our DE results are driven by labor rather than covid status.

The TIL vs TNL differential gene expression results from (Pique-Regi et al. 2019) can be used for this. I downloaded the supplementary file. The file is described as follows: “The columns represent: 1) Cluster or cell-type name, 2) Comparison groups or contrast (i.e., TNL vs TIL, TIL vs PTL), 3) Ensembl gene identifier, 4) Gene symbol, 5) baseMean gene baseline expression as calculated by DESeq2, 6) log2 Fold Change of the first group in column two versus the second group, 7) Standard error estimated for the log2 Fold Change, 8) Nominal p-value, 9) q-value or adjusted p-value to control for FDR. Only rows with q < 0.2 are reported.”.

til <- readxl::read_excel("../info/from-papers/pique-regi_2019/supp/elife-52004-supp5-v2.xls")
labor <- til %>% filter(Compare == "TNL_v_TIL" & padj < 0.05)
head(labor)

Now we have to find out how many of the DE genes between TNL and TIL are shared with DE genes between covid and control.

# all DE genes across all cell types between covid and control
de.all <- de.genes

for(i in names(de.all)) {
  de.all[[i]]$gene <- rownames(de.all[[i]])
}

de.all <- do.call(rbind, de.all)

de.all.genes <- de.all$gene %>% unique()

# all labor DE genes
labor.genes <- labor$Gene %>% unique()

# number of covid vs control DE genes (across all cell types)
length(de.all.genes)

[1] 1267

# number of TNL vs TIL DE genes (across all cell types)
length(labor.genes)

[1] 195

# number of DE genes shared between covid vs control and TNL vs TIL
length(intersect(de.all.genes, labor.genes))

[1] 84

There are 84 genes that are shared between covid vs control and TNL vs TIL. These are likely related to labor. These genes are as follows:

intersect(de.all.genes, labor.genes) %>% sort()

 [1] "ABCA1"    "ACTB"     "AHR"      "ANXA1"    "B4GALT1"  "C1QA"     "C1QB"    
 [8] "C1QC"     "CCL18"    "CCND2"    "CCNL1"    "CD163"    "CD2"      "CD248"   
[15] "CD3E"     "CD55"     "CDKN2A"   "CFLAR"    "CKLF"     "CLDN5"    "CXCL1"   
[22] "DEFB1"    "DNAJA4"   "ELF1"     "ELF3"     "EMP1"     "ETS2"     "FABP5"   
[29] "FCGR2B"   "FILIP1L"  "FNBP1"    "FNIP2"    "FOSB"     "GNLY"     "GPCPD1"  
[36] "HLA-DRB5" "HSD11B1"  "HSPA6"    "ICAM1"    "ID2"      "IDH2"     "IFI27"   
[43] "IFITM2"   "IFITM3"   "IGFBP1"   "IGFBP2"   "IGFBP4"   "IGFBP5"   "IL1RL1"  
[50] "JUND"     "KCNQ1OT1" "KDM6B"    "KLF6"     "LMNA"     "LSP1"     "MAFF"    
[57] "MALAT1"   "MBNL1"    "MGST1"    "NAMPT"    "NEAT1"    "NFKB1"    "NFKBIA"  
[64] "PDE4B"    "PNRC1"    "POU2F2"   "PTGDS"    "RASGEF1B" "RGCC"     "RPS26"   
[71] "SAMSN1"   "SERPINE1" "SOD2"     "SPP1"     "SRGN"     "SRSF5"    "STAB1"   
[78] "TFPI"     "TIPARP"   "TNFAIP2"  "TNFAIP3"  "TXN"      "TYMS"     "ZFP36"   

Now we can compare whether the shared genes change in the same direction that we expect based on the information that our control samples come from C-sections.

# create a dataframe of shared genes and their direction of changes
shared <- data.frame(
  gene = intersect(de.all.genes, labor.genes) %>% sort()
)

shared$labor_logfc <- NA
for(i in 1:nrow(shared)){
  shared$labor_logfc[i] <- mean(labor$log2FoldChange[labor$Gene == shared$gene[i]])
} # because the same gene can be DE in multiple cell types, we have to take average logFC. Not the best option, but quick and dirty. 

shared$covid_logfc <- NA
for(i in 1:nrow(shared)){
  shared$covid_logfc[i] <- mean(de.all$avg_logFC[de.all$gene == shared$gene[i]])
} 

shared$labor_pattern[shared$labor_logfc < 0] <- "TIL_high"
shared$labor_pattern[shared$labor_logfc > 0] <- "TNL_high"

shared$covid_pattern[shared$covid_logfc < 0] <- "control_high"
shared$covid_pattern[shared$covid_logfc > 0] <- "covid_high"

The contigency table of shared DE genes looks like this:

table(shared$labor_pattern, shared$covid_pattern)

          
           control_high covid_high
  TIL_high            9         59
  TNL_high            3         13

Or in terms of proportions, like this:

table(shared$labor_pattern, shared$covid_pattern) %>% prop.table()

          
           control_high covid_high
  TIL_high   0.10714286 0.70238095
  TNL_high   0.03571429 0.15476190

Fisher’s exact test for this contingency table looks like this:

table(shared$labor_pattern, shared$covid_pattern) %>% fisher.test(alternative = "two.sided")

    Fisher's Exact Test for Count Data

data:  .
p-value = 0.6912
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.1383802 4.3366594
sample estimates:
odds ratio 
  0.664588 

Labor is likely a confounding factor in our analysis, but the effect is very mild, if any. Since our control samples are from C-sections (i.e. not in labor), but covid samples are at term and in labor (not C-section), if labor is a confounding factor, we should expect a high proportion of genes to be shared between “covid_high” and “TIL_high”. We see this looking at the odds ratio only. But the p value for Fisher’s exact test is 0.69, suggesting that the labor and covid status are largely independent in our samples.

While it’s nice to be aware of this potential but mild effect, it is not a big problem. This effect is only relevant for genes that are shared between the two comparisons—TNL vs TIL and covid vs control. That fraction of genes, to begin with, is small. Around 195 genes are DE between TIL and TNL. Of these 84 are also DE between covid and control, likely representing the “labor signature” rather than “covid signature”. But the number of genes that are DE between covid and control is much larger (1267), i.e. over 1000 genes are DE between covid and control that are not DE between TNL and TIL. These most likely represent the “covid signature” that we are interested in.

8 Data-vis ideas for the paper

8.1 DE genes dotplot

We need to find a good way to present DE genes. Violin plots are nice if the number of to show is small. If we want to show a small set of genes, we can use violin plots.

Another option is dotplots. The advantage of dotplots is that they are much more space-efficient, which is desirable for paper figures. In a 6inx8in figure, we can comfortably show data for about 50 genes in all cell types separated by covid status.

The Seurat::DotPlot function has some problems, though, that make is not really usable for good quality plotting. One big problem is how it handles splitting samples by treatment and control (covid and control in our case) in split.by argument. It does split your cell types by treatment status so that the average expression by treatment status for each cluster can be visualized, but it uses two different colours for control and treatment. This causes two problems. First, ggplot doesn’t support multiple colour scales, and Seurat, instead of finding a way around the issue, has chosen to not show the colourbar guide at all when split.by is used (for more on this and another related bug, see https://github.com/satijalab/seurat/issues/2487 and https://github.com/satijalab/seurat/issues/2342.). Second, using two colors for values (avg.exp) on the same scale is a bad design choice, especially when the point is to compare them side by side—intensities for two colours mapping to the exact same avg.exp value can be perceived differently just because they are two different colors.

1. We have to do a work-around to plot the dotplots split by covid status, plotting with a single color. We just use ggplot directly to make the dotplots, with covid status on one axis, genes on another, and faceted by cell type. This requires us to do the calculations for avg.exp.scaled and pct.exp ourselves. For this we can cheat—we can run Seurat::dotplot with celltype_covid as idents so it calculate those values split by covid status for each cell type, but instead of plotting it, we can just extract the $data slot from the resulting ggplot object and use it for our custom plotting.
2. We can also show result of DE analysis on top of this plot. We can add a line segment connecting the covid and control dots in those cell types in which the p value of DE test of below threshold.
3. For gene order, we can use hclust for clustering, and use that order for sorting genes. This way, genes that have similar expression patterns are positioned close to each other, and there is some structure to the plotting order.

Here for demo, we’ll use top 50 DE genes from dec.APC.

# theme for DE genes faceted dotplot
theme_dotplot <- theme(
  panel.background = element_blank(),
  text = element_text(color = "Black"),
  panel.grid = element_blank(),
  panel.border = element_rect(fill = NA, size = 0.1, color = "grey60"),
  axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6, color = "Black"),
  axis.text.y = element_text(size = 6, color = "Black"),
  axis.title = element_blank(),
  axis.ticks = element_line(size = 0.1, color = "Black"),
  strip.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5, size = 6.5, color = "Black"),
  strip.background = element_blank(),
  legend.title = element_text(size = 6),
  legend.position = "right",
  legend.box.spacing = unit(0.1, "lines"),
  legend.background = element_blank(),
  legend.title.align = 0,
  legend.text = element_text(size = 5.5),
  legend.key.size = unit(0.5, "lines"),
  legend.key = element_blank(),
  legend.spacing.y = unit(0.1, "lines"),
  panel.spacing.x = unit(0, "lines")
)

# extract plot data from Seurat::DotPlot
Idents(seur) <- seur$celltype_covid
p <- Seurat::DotPlot(object = seur, assay = "SCT", 
                     features = rownames(de.genes$dec.APC)[1:50])

p.dat <- p$data

p.dat$celltype <- gsub(p.dat$id, pattern = "(.+)_([a-z]{5})", replacement = "\\1")
p.dat$covid <- gsub(p.dat$id, pattern = "(.+)_([a-z]{5})", replacement = "\\2")

# clustering for ordering
cdat <- pivot_wider(data =  p.dat, id_cols = 3, names_from = id, values_from = avg.exp.scaled)

cdat <- as.data.frame(cdat)
rownames(cdat) <- cdat$features.plot
cdat$features.plot <- NULL

cor.matrix <- cor(t(cdat))

# reorder correlation matrix based on clustering
dd <- as.dist((1 - cor.matrix)) 
hc <- hclust(dd, method = 'complete')
cdat  <- cdat[hc$order, ]

p.dat$features.plot <- factor(p.dat$features.plot, levels = rownames(cdat))

# prepare annotation data to display DE stats. This gives a df of only the combination of celltypes and genes where adj.p < 0.05

de.ann <- de.genes

for(i in names(de.ann)) {
  de.ann[[i]]$gene <- rownames(de.ann[[i]])
  de.ann[[i]]$celltype <- i
} # add gene names as column. the celltype info too, so that we can rbind all list items.

de.ann <- do.call(rbind, de.ann)
de.ann <- subset(de.ann, gene %in% p.dat$features.plot)
de.ann$diff.exp <- de.ann$p_val_adj < 0.05
de.ann$diff.exp <- tolower(as.character(de.ann$diff.exp))

de.ann$p_val <- NULL
de.ann$avg_logFC <- NULL
de.ann$pct.covid <- NULL
de.ann$pct.cntrl <- NULL
rownames(de.ann) <- NULL

q <- ggplot(data = p.dat, aes(x = covid, y = features.plot)) +
  geom_point(aes(fill = avg.exp.scaled, size = pct.exp), shape = 21, stroke = 0, 
             color = "White") +
  scale_fill_gradient(low = "White", high = "Red", name = "avg.exp\n(scaled)") +
  scale_x_discrete(breaks = c("cntrl", "covid"), 
                   labels = c("Control", "Covid")) +
  scale_size_area(max_size = 3.5,
                  guide = guide_legend(
                    override.aes = list(
                      color = "White", 
                      fill = "Black"))
                  ) +
  facet_grid(. ~ celltype) +
  ungeviz::geom_hpline(data = de.ann,
                       aes(y = gene, x = 1.5, color = diff.exp),
                       size = 0.35, width = 1, lineend = "butt") +
  scale_color_manual(breaks = "true",
                     labels = "< 0.05",
                     values = "Black",
                     name = "DE adj.p") +
  theme_dotplot


cowplot::ggsave2(q, 
                 filename = "../results/04_de-genes-by-celltype/logfc_0.40/de/plots/dotplot_top50-dec.APC.pdf", 
                 width = 6, height = 7, units = "in")

print(q)

8.2 GO enrichment dotplot

We also have to find a way to visualize the GO enrichment results. This one is tricky, but dotplot of enrichment and p value is a good option. But there are too many GO categories that enriched across all celltypes (around 1000), so we do not to filter them somehow.

8.2.1 Calculate enrichment

First we put all GO enrichment results in the same dataframe and then calculate the enrichment. Here the GO categories will be repeated because the same GO category can be enriched in multiple cell types, but they will of course have different enrichment and p values in each cell type.

# prepare data
go.pdat <- de.go[lapply(de.go, nrow) > 0]

for(i in names(go.pdat)) {
  go.pdat[[i]]$celltype <- i
}

go.pdat <- do.call(rbind, go.pdat)

go.pdat$pvalue <- NULL
go.pdat$geneID <- NULL
rownames(go.pdat) <- NULL

# calculate enrichment
go.pdat$b <- as.numeric(sapply(
    strsplit(go.pdat$GeneRatio, split = "/"), MARGIN = 1, FUN = '[['))

go.pdat$n <- as.numeric(sapply(
    strsplit(go.pdat$GeneRatio, split = "/"), MARGIN = 2, FUN = '[['))

go.pdat$B <- as.numeric(sapply(
    strsplit(go.pdat$BgRatio, split = "/"), MARGIN = 1, FUN = '[['))

go.pdat$N <- as.numeric(sapply(
    strsplit(go.pdat$BgRatio, split = "/"), MARGIN = 2, FUN = '[['))

go.pdat$enrichment <- (go.pdat$b / go.pdat$n) / (go.pdat$B / go.pdat$N)

head(go.pdat)

8.2.2 Plotting function

Since the same GO categories can be enriched in different celltypes, to avoid repetition and to enable seeing patterns of enrichment across celltypes, we want to visualize enrichment output not individually by celltype, but in a combined plot for all celltypes. We can dotplot celltype on one axis, and GO category on another. Color of the dot represents p value and size represents enrichment. We can cluster the GO categories by and order them accordingly so that we can easily see the patterns across celltypes.

# theme for GO dotplot
theme_GO <- theme(
  panel.background = element_blank(),
  text = element_text(color = "Black"),
  panel.grid.major = element_line(size = 0.15, color = "Grey90"),
  panel.border = element_blank(),
  axis.line = element_line(size = 0.25, color = "Black"),
  axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6, color = "Black"),
  axis.text.y = element_text(size = 6, color = "Black"),
  axis.title = element_blank(),
  axis.ticks = element_line(size = 0.1, color = "Black"),
  strip.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5, size = 6.5, color = "Black"),
  legend.title = element_text(size = 6),
  legend.direction = "vertical",
  legend.box = "vertical",
  legend.box.spacing = unit(0.2, "lines"),
  legend.box.just = "left",
  legend.position = "right",
  legend.background = element_blank(),
  legend.title.align = 0,
  legend.text = element_text(size = 5.5),
  legend.key.size = unit(0.5, "lines"),
  legend.key = element_blank(),
  legend.spacing.y = unit(0.2, "lines"),
)

# function

dotplot_go <- function(go.pdat.sub) {
  
  # clustering for ordering
  cdat <- pivot_wider(data = go.pdat.sub, id_cols = 2, 
                      names_from = celltype, values_from = p.adjust)
  cdat[, 2:ncol(cdat)] <- apply(cdat[, 2:ncol(cdat)], 
                                MARGIN = 2, FUN = function(x){replace_na(x, 1)})
  
  cdat <- as.data.frame(cdat)
  rownames(cdat) <- cdat$Description
  cdat$Description <- NULL
  
  cor.matrix <- cor(t(cdat))
  
  dd <- as.dist((1 - cor.matrix)) 
  hc <- hclust(dd, method = 'complete')
  
  cdat <- cdat[hc$order, ]
  
  go.pdat.sub$Description <- factor(go.pdat.sub$Description, 
                                    levels = rownames(cdat))
  # plot
  p <- ggplot(data = go.pdat.sub,
              aes(x = celltype, y = Description)) +
    geom_point(aes(size = enrichment, color = -log10(p.adjust))) +
    scale_color_gradient(low = colorRampPalette(c("White", "Blue"))(10)[5],
                         high = colorRampPalette(c("White", "Blue"))(10)[10]) +
    scale_size_area(max_size = 3.5) +
    theme_GO
  
  return(p)
}

8.2.3 Dotplot without redundancy filtering

Here, after having combined all enriched GO categories from all cell types, we sort them by p value and pick the top 100, regardless of which cell type they come from. This is perhaps not the best way to filter, but one option that picks up the most highly enriched categories across most celltypes. But it does miss some celltypes whose enriched GO categories are not in the top 100.

# plotting subset
go.sub <- go.pdat[order(go.pdat$p.adjust, decreasing = FALSE)[1:100], ]

p <- dotplot_go(go.pdat.sub = go.sub)

ggsave(
  p, 
  filename = "../results/04_de-genes-by-celltype/logfc_0.40/go/plots/go_dotplot_top100.pdf",
  width = 5, height = 5, units = "in"
)

print(p)

8.2.4 Plotting after REVIGO

REViGO (http://revigo.irb.hr/) is an online tool that takes a list of GO terms and p values as input and outputs a reduced list of GO categories. It does this by finding GO terms that are similar to each other, and based on shared genes, it assigns “uniqueness” and “dispensability” scores for each GO category. We can filter out GO categoeis that have high dispensability score—these GO categories are redundant because they have been represented by other GO categories in the list.

# write for running revigo
write.table(
  go.pdat[, c("ID", "p.adjust")], 
  file = "../results/04_de-genes-by-celltype/logfc_0.40/go/revigo/GO_for-revigo.txt", 
  quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)

revigo <- read.csv("../results/04_de-genes-by-celltype/logfc_0.40/go/revigo/REVIGO.csv")

# plotting subset
to.keep <- revigo$term_ID[revigo$dispensability < 0.85]
go.sub <- go.pdat[go.pdat$ID %in% to.keep, ]

p <- dotplot_go(go.pdat.sub = go.sub)

ggsave(
  p, 
  filename = "../results/04_de-genes-by-celltype/logfc_0.40/go/plots/go_dotplot_revigo.pdf",
  width = 5.75, height = 6, units = "in"
)

print(p)

8.2.5 Plotting after REViGO: only terminal nodes of DAG

GO terms are hierarchically organized in a directed acyclic graph (DAG). Sometimes you can have two terms enriched where one is more inclusive than the other, and the less inclusive one is the more informative one. For example, if you have “cell killing” and “negative regulation of cell killing” as two enriched GO terms, the latter—the more informative one—is a child term of the the former. The parent-child relationship of GO terms is another axis along which we can prune the list of GO terms, by discarding terms whose child terms are also enriched.

# function for extracting terminal nodes (terms whose child terms are not enriched)
terminal <- function(terms, ontology=c("C", "P", "F"))
{
  FUN <- switch(match.arg(ontology),  
                C = GOCCPARENTS,
                P = GOBPPARENTS, 
                F = GOMFPARENTS)
  terminal <- terms
  seen <- c(terms, "all")
  while (length(terms)) {
    seen <- c(terms, seen)
    terms <- mappedRkeys(FUN[terms])
    terminal <- terminal[!terminal %in% terms]
    terms <- terms[!terms %in% seen]
  }
  terminal
}

# run revigo output terms through terminal function
leaf <- terminal(terms = as.character(revigo$term_ID), ontology = "P")
length(leaf)

[1] 46

# plot
p <- dotplot_go(go.pdat.sub = go.pdat[go.pdat$ID %in% leaf, ])

ggsave(
  p, 
  filename = "../results/04_de-genes-by-celltype/logfc_0.40/go/plots/go_dotplot_revigo_terminal-nodes.pdf",
  width = 6, height = 5.25, units = "in"
)

print(p)

9 Session Info

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] GO.db_3.11.4           org.Hs.eg.db_3.11.4    AnnotationDbi_1.50.3  
 [4] IRanges_2.22.2         S4Vectors_0.26.1       Biobase_2.48.0        
 [7] BiocGenerics_0.34.0    clusterProfiler_3.16.1 ungeviz_0.1.0         
[10] ggtext_0.1.0           cowplot_1.1.0          ggforce_0.3.2         
[13] ggrepel_0.8.2          ggthemes_4.2.0         patchwork_1.0.1       
[16] Seurat_3.2.1           forcats_0.5.0          stringr_1.4.0         
[19] dplyr_1.0.2            purrr_0.3.4            readr_1.3.1           
[22] tidyr_1.1.2            tibble_3.0.3           ggplot2_3.3.2         
[25] tidyverse_1.3.0       

loaded via a namespace (and not attached):
  [1] backports_1.1.10      readxl_1.3.1          fastmatch_1.1-0      
  [4] plyr_1.8.6            igraph_1.2.5          lazyeval_0.2.2       
  [7] splines_4.0.2         BiocParallel_1.22.0   listenv_0.8.0        
 [10] urltools_1.7.3        digest_0.6.25         htmltools_0.5.0      
 [13] GOSemSim_2.14.1       viridis_0.5.1         fansi_0.4.1          
 [16] magrittr_1.5          memoise_1.1.0         tensor_1.5           
 [19] cluster_2.1.0         ROCR_1.0-11           limma_3.44.3         
 [22] globals_0.12.5        graphlayouts_0.7.0    modelr_0.1.8         
 [25] matrixStats_0.56.0    enrichplot_1.8.1      prettyunits_1.1.1    
 [28] colorspace_1.4-1      rvest_0.3.6           blob_1.2.1           
 [31] haven_2.3.1           xfun_0.16             crayon_1.3.4         
 [34] jsonlite_1.7.1        scatterpie_0.1.4      spatstat_1.64-1      
 [37] spatstat.data_1.4-3   survival_3.2-3        zoo_1.8-8            
 [40] glue_1.4.2            polyclip_1.10-0       gtable_0.3.0         
 [43] leiden_0.3.3          future.apply_1.6.0    abind_1.4-5          
 [46] scales_1.1.1          DOSE_3.14.0           DBI_1.1.0            
 [49] miniUI_0.1.1.1        Rcpp_1.0.5            gridtext_0.1.1       
 [52] viridisLite_0.3.0     xtable_1.8-4          progress_1.2.2       
 [55] gridGraphics_0.5-0    reticulate_1.16       bit_4.0.4            
 [58] europepmc_0.4         rsvd_1.0.3            htmlwidgets_1.5.1    
 [61] httr_1.4.2            fgsea_1.14.0          RColorBrewer_1.1-2   
 [64] ellipsis_0.3.1        ica_1.0-2             pkgconfig_2.0.3      
 [67] farver_2.0.3          dbplyr_1.4.4          uwot_0.1.8           
 [70] deldir_0.1-28         labeling_0.3          ggplotify_0.0.5      
 [73] tidyselect_1.1.0      rlang_0.4.7           reshape2_1.4.4       
 [76] later_1.1.0.1         munsell_0.5.0         cellranger_1.1.0     
 [79] tools_4.0.2           cli_2.0.2             downloader_0.4       
 [82] generics_0.0.2        RSQLite_2.2.0         broom_0.7.0          
 [85] ggridges_0.5.2        fastmap_1.0.1         yaml_2.2.1           
 [88] goftest_1.2-2         fs_1.5.0              knitr_1.30           
 [91] bit64_4.0.4           fitdistrplus_1.1-1    tidygraph_1.2.0      
 [94] RANN_2.6.1            ggraph_2.0.3          pbapply_1.4-3        
 [97] future_1.18.0         nlme_3.1-149          mime_0.9             
[100] DO.db_2.9             xml2_1.3.2            compiler_4.0.2       
[103] rstudioapi_0.11       plotly_4.9.2.1        png_0.1-7            
[106] spatstat.utils_1.17-0 reprex_0.3.0          tweenr_1.0.1         
[109] stringi_1.5.3         strapgod_0.0.4.9000   lattice_0.20-41      
[112] Matrix_1.2-18         vctrs_0.3.4           pillar_1.4.6         
[115] lifecycle_0.2.0       BiocManager_1.30.10   triebeard_0.3.0      
[118] lmtest_0.9-38         RcppAnnoy_0.0.16      data.table_1.13.0    
[121] irlba_2.3.3           httpuv_1.5.4          qvalue_2.20.0        
[124] R6_2.4.1              promises_1.1.1        KernSmooth_2.23-17   
[127] gridExtra_2.3         codetools_0.2-16      assertthat_0.2.1     
[130] MASS_7.3-52           withr_2.3.0           sctransform_0.2.1    
[133] mgcv_1.8-31           hms_0.5.3             grid_4.0.2           
[136] rpart_4.1-15          rvcheck_0.1.8         Rtsne_0.15           
[139] lubridate_1.7.9       shiny_1.5.0          

10 References

Pique-Regi, R., R. Romero, A. L. Tarca, E. D. Sendler, Y. Xu, V. Garcia-Flores, Y. Leng, F. Luca, S. S. Hassan, and N. Gomez-Lopez. 2019. “Single Cell Transcriptional Signatures of the Human Placenta in Term and Preterm Parturition.” eLife 8: e52004.

Rusinova, I., S. Forster, S. Yu, A. Kannan, M. Masse, H. Cumming, R. Chapman, and P. J. Hertzog. 2012. “INTERFEROME v2.0: an updated database of annotated interferon-regulated genes.” Nucleic Acids Research 41 (D1): D1040–D1046. https://doi.org/10.1093/nar/gks1215.